
from pylab import *
import scipy.misc
import scipy.ndimage

ion()
# aa = array([[-.0147, -.0933, -.1253, -.0933, -.0147],
#             [-.0469, -.0640, -.0640, -.0640, -.0469],
#             [0,0,0,0,0],
#             [ .0469,  .0640,  .0640,  .0640,  .0469],
#             [ .0147,  .0933,  .1253,  .0933,  .0147]]).T / 1.9368

aa = array([[-.0147, -.0469, 0, .0469, .0147],
            [-.0933, -.0640, 0, .0640, .0933],
            [-.1253, -.0640, 0, .0640, .1253],
            [-.0933, -.0640, 0, .0640, .0933],
            [-.0147, -.0469, 0, .0469, .0147]]) / 1.9368

so = array([[-1.,-2,-1],
            [0,0,0],
            [1,2,1]])/8.

bb=scipy.array([
      [ -0.003776, -0.010199, 0., 0.010199, 0.003776 ],
      [ -0.026786, -0.070844, 0., 0.070844, 0.026786 ],
      [ -0.046548, -0.122572, 0., 0.122572, 0.046548 ],
      [ -0.026786, -0.070844, 0., 0.070844, 0.026786 ],
      [ -0.003776, -0.010199, 0., 0.010199, 0.003776 ]
      ])

qq = scipy.misc.lena()


dxa = scipy.ndimage.convolve(qq, aa)
dya = scipy.ndimage.convolve(qq, aa.T)
dxb = scipy.ndimage.convolve(qq, bb)
dyb = scipy.ndimage.convolve(qq, bb.T)

mma = sqrt(dxa**2+dya**2)
mmb = sqrt(dxb**2+dyb**2)


vv=40
figure(1)
subplot(2,3,1)
imshow(qq, cmap=cm.bone, vmin=0, vmax=255, interpolation='nearest')
subplot(2,3,2)
imshow(dxa, cmap=cm.RdBu, vmin=-vv, vmax=vv, interpolation='nearest')
subplot(2,3,3)
imshow(dya, cmap=cm.RdBu, vmin=-vv, vmax=vv, interpolation='nearest')

subplot(2,3,5)
imshow(dxb, cmap=cm.RdBu, vmin=-vv, vmax=vv, interpolation='nearest')
subplot(2,3,6)
imshow(dyb, cmap=cm.RdBu, vmin=-vv, vmax=vv, interpolation='nearest')

figure(2)
vv = .2
subplot(1,2,1)
title('Zernike A11')
imshow(aa, cmap=cm.RdBu, vmin=-vv, vmax=vv, interpolation='nearest')
subplot(1,2,2)
title('Shigeru')
imshow(bb, cmap=cm.RdBu, vmin=-vv, vmax=vv, interpolation='nearest')

figure(3)
subplot(2,2,1)
imshow(qq, cmap=cm.bone, vmin=0, vmax=255, interpolation='nearest')
axis([250,500,400,150])
subplot(2,2,3)
imshow(mma*(arctan2(dya,dxa)+pi+.5), cmap=cm.RdBu, vmin=-9, vmax=9, interpolation='nearest')
axis([250,500,400,150])
subplot(2,2,4)
imshow(mmb*(arctan2(dyb,dxb)+pi+.5), cmap=cm.RdBu, vmin=-9, vmax=9, interpolation='nearest')
axis([250,500,400,150])
